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Abstract 

In this paper we exploit a novel idea for the optimization of flows governed by the 
Euler equations. The algorithm consists of marching on the design hypersurface while 
improving the distance to the state and costate hypersurfaces. We consider the problem 
of matching the pressure distribution to a desired one, subject to the Euler equations, 
both for subsonic and supersonic flows. The rate of convergence to the minimum for 
the cases considered is 3 to 4 times slower than that of the analysis problem. Results 
are given for Ringleb flow and a shockless recompression case. 
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1 Introduction 


In recent years there has been a renewed interest in optimal design in fluid mechanics. 
Faster computers and reliable numerical simulations make feasible some of the aerodynamics 
optimal design problems which are of engineering interest. 

The statement above becomes only partially true when we consider either flows governed 
by Euler or Navier-Stokes equations, or complicated geometrical configurations with many 
control parameters. For such cases, shape optimization seems to be still not practical due 
to extremely time consuming computation. The present work aims at a flexible and feasible 
approach for such intensive computational problems by applying a novel algorithm proposed 
by Ta’asan in [19]. 

The problem of finding a shape that achieves given performance has been attacked by 
means of inverse problem formulations [12], [2], [21], [5]. These methods have in common the 
advantage of being solved at the same cost of an analysis problem. They are in general not 
extendable to three dimensions. Moreover the set of problems that can be solved by means 
of inverse design is limited. 

A more general framework is to consider aerodynamics design problems as optimization 
problems. From the mathematical viewpoint the problem is to find U such that 

f UeU 

\ £{U) < £(V) W g U 

where U is a given set and £ is a real-valued functional defined on U . 

Shape design optimization problems are tightly related to control of a system governed 
by partial differential equations where the controls are on the boundary. Lions set in [13] 
the mathematical framework for such problems. The theory is concerned mainly with linear 
systems and is devoted “(i) to obtain necessary (or possibly necessary and sufficient) con- 
ditions for U to be an extremum (or minimum), (ii) to study the structure and properties of 
the equations expressing these conditions, (iii) to obtain constructive algorithms amenable 
to numerical computations for the approximation of U v . 

Pironneau ([15], [16]) derived an adjoint method for the minimum-drag problem in Stokes 
flows and subsequently in flows governed by the incompressible Navier-Stokes equations. 
Since a Navier-Stokes solver was not available, some solutions were obtained using simpler 
models; see Glowinski and Pironneau [6]. 

For the Euler equations Jameson proposed in [10] an adjoint method for wing design 
which makes use of conformal mapping to control the shape of the wing. Iollo, Salas and 
Ta’asan [9] studied the case of Euler flows with embedded shocks for a one-dimensional 
case, and discussed the boundary conditions for the adjoint equations. At the shock location 
it was shown that further conditions are needed for the adjoint equation to be well posed. 
Subsequently, Iollo and Salas extended these results to two-dimensional flows, and presented 
computations with higher-order spatial accuracy [8]. 

The high computational cost for solving optimization problems governed by fluid dynamics 
equations comes from several sources. The first is the cost of a single analysis, which for the 
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Navier-Stokes equations in three dimensions is of the order of a few CRAY hours. Another 
source is the fact that a repeated solution of the flow equations may be required for gradient 
methods. An additional significant cost may arise from the calculation of the gradient of 
the functional. 

The use of adjoint methods eliminates the unnecessary cost resulting from computation 
of the gradient, and is much more efficient compared to other methods including finite- 
differences and sensitivity analysis. It requires the computation of an extra system of 
partial differential equations (PDEs), namely the costate equations, but the total cost for 
gradient calculation is independent of the number of design variables. A comparison study 
of calculating gradients using adjoint methods and Unite-difference methods was done by 
Beux and Dervieux [3]. They also solved pressure reconstruction problems for compressible 
internal flows, comparing the performances of several algorithms. Flows with embedded 
shocks were not considered in this work. 

The adjoint method, being an efficient method for calculating the gradient, does not ad- 
dress the computational expense related to the number of gradient iterations required to 
reach the minimum. In general, the number of iterations required to achieve the minimum 
grows more than linearly with the number of controls used, making infeasible design prob- 
lems in three dimensions with many design variables. 

Ta’asan proposed in [18] an algorithm to reduce the cost of the optimization to that of a 
single analysis, namely the one shot method. The idea is to solve the flow equations, the 
costate equations and the optimality condition at the same time. The main idea in that 
algorithm was to perform the optimization iteration on coarse grids that are used anyway 
in the multigrid process. Small numbers of design variables were considered in that case. 

Ta’asan, Kuruvila and Salas [20] applied this technique to a potential flow, and extended 
the method to cases of moderate numbers of design variables. Different design variables 
are associated with different grids depending on the smoothness of the shape functions 
associated with them, and are updated on these grids. The performance of this algorithm 
was practically independent of the number of design variables. 

Arian and Ta’asan [1] extended the one shot method to infinite-dimensional design space. 
The main idea of the method was to construct a relaxation that smoothes the errors in 
the design variables. Application to control problems and shape design problems have 
demonstrated solution of the full optimization problem in a cost comparable to that of solving 
the analysis just a few times, independent of the number of design variables (experiments 
using up to 128 design variables have been done). 

Beux and Dervieux [4] proposed a hierarchical strategy in which the number of control 
parameters is progressively increased performing a multilevel optimization that seems to 
render the computational cost independent from the number of control parameters. 

The drawbacks of the one shot methods are their programming complexity and the fact 
that their use is limited to multigrid solvers. This was the motivation for the study of a new 
type of solution strategy for optimization problems governed by PDEs [19]. The goal was to 
try to get methods that solve the optimization problem in a cost comparable to that of the 
analysis. The emphasize was on simplicity and flexibility to work in existing frameworks 
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which do not necessarily involve multigrid methods. 

The main observation is the following. The solution of the optimization problem lies on 
the intersection of the state, costate and design hypersurfaces (in the state, costate and 
design spaces). Gradient-based methods (including adjoint formulations) can be viewed as 
marching along the intersection of the state and costate hypersurfaces. This is an expensive 
process since each step requires the solution of two PDEs. The idea of the pseudo-time 
method was to perform the marching on the design hypersurface while improving the dis- 
tance to the other two. The cost of such an iteration per step is significantly smaller than 
that of gradient-based methods. Its convergence has been shown by Ta’asan in [19] to be 
independent of the number of design variables. 

In the present paper we apply the pseudo-time method to optimization problem using the 
Euler equations. Using this method the cost of optimization becomes of the same order as 
that of analysis. Moreover the algorithm may be implemented with no substantial changes 
to existing codes. Numerical results indicate that the method converges at a rate which is 
independent of the number of design variables. 


2 Problem statement 

The Euler equations are given by 


where 


with 


Uf + + G y — 0 
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p = density 

u = x-component of the velocity vector 
v = y-component of the velocity vector 
e = specific total energy 
p = pressure 
a = speed of sound 
7 = ratio of specific heats 

7 - 1 


and p = up(2e — u 2 — v 2 ). Furthermore let 
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We assume that these equations are defined on a domain $ which includes a sub-domain 
0 whose boundary is denoted by L. On the boundary we define a curvilinear coordinate s 
and a normal n = (n X} n y ) pointing outward. A real valued functional £(T, V (T)) is given, 
where V(T) is the solution of the Euler equations with boundary conditions on L. The 
optimization problem that we study is 

minimize the functional £(T, V(T)) over all the admissible shapes of the boundary L. 

We focus on the following model problem. The sub-domain 0 is represented by a nozzle; 
see Fig. 1. At the inlet, total pressure, total temperature and the ratio a = v/u are fixed. 
At the outlet, if the flow is subsonic, the static pressure is fixed and at the solid walls the 
impermeability condition un x + vn y = 0 is enforced. The upper wall is kept fixed. The 
lower wall 0 is represented by mean of the parameterization 

y(@) = ( 3 ) 

i 


where the functions fi(x) are some shape functions and a = (aq, ..., cq, ...) is the corres- 
ponding set of shape coefficients. Given a desirable lower wall pressure distribution p*(x) 
and denoting by p e (x) the actual one on the lower wall, the optimization problem consists 
in finding a set of shape coefficients cq such that the functional 

£ = 1, J a (p & ~ P*) 2dx (4) 


is minimized. 


3 Optimality conditions 


The optimality conditions are derived by introducing Lagrange multipliers and considering 
the augmented functional 

£(U,a,A,/i) = £(a,U) + f *A(AU* + BUjdfi + f /ip\ ■ nds (5) 

Jn Je 

where V = (it, v). The vector A(x, y) = * ( Ai , A 2 , A 3 , A 4 ), and the scalar p(s) are the Lagrange 
multipliers. 

Calculating the variation of the functional C with respect to the variation of the functions 
U, A, fi and the parameters cq respectively, we obtain (see [8]) 


[ b dp 

SCu = i au 


( p & — p*)XJdx + / f A(An x + Hn y )\Jds + 
j r 


A + 'A,B)UdO + J & pn ^ 


U ds 


( 6 ) 
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where 


dp -2k ( U2 + V2 -u -v l) and ^ 0 1 0 0 

— Zfh l _ ? — i — V ? 1 | dilU. 
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( 7 ) 

( 8 ) 


(9) 


where 0 is the angle between the normal n and the y-axis, t = (—n y , n x ) } and U, A, p and 
cq are the variations of the corresponding arguments. 

At the minimum of the functional, for all the possible choices of the functions U, A, p 
and of the parameters a, we must have 

SCu = SCa = SC^ = SC a = 0. (10) 


Therefore, we have 


SjCa = 0 AU-c + BUj = 0 on 0 


and 

SCfj, = 0 pV • n = 0 on 0 


which are the Euler equations and boundary conditions. Furthermore 

SCu = 0 f AA x + t BA y = 0 on 0 


( 11 ) 


and 


where 


dp 

dU 


(p 


0 


p*) cos 9 + t A(A n x + B n y ) + pn 


dpV 

~dV 


0 on 0 


( 12 ) 


f. i — — Ai -|- u\2 — kV 2 )X 4 


( 13 ) 


For the boundary condition on inlet and outlet we refer the reader to [8]. Given U, the set 
of costate eqs. (11-13) determine uniquely A in 0 and p on 0. Finally, given a and knowing 
U and A, we can calculate from eq.(9) 
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+ f /a — ^ — - • n fids — f /ipY ■ t cos 2 0 ds + 

J © oy Je ax 

+ f ypV • n —j-— sin 0 cos 6 ds (14) 

Je ax 

In case of a shock occuring in the flow held, we split the domain of integration by means of 
a curve T that coincides with the shock where it exists. Then we follow the same derivation 
presented so far on each of the two sub-domains, regarding T as a boundary; see [8]. The 
resulting extra condition for A on the shock is A = 0. It should be noted that if the 
shocks are not treated properly, the problem of solving the costate equations with boundary 
conditions is not well-posed. Jameson presented in [11] results for transonic hows over 
airfoils where the wave-drag is minimized. He does not use any special treatment for the 
shock but the costate equations converge. This is due to the fact that the scheme that 
he uses for solving the Euler equations smears the shocks over several grid points, due to 
artihcial viscosity. 


4 Pseudo-time optimization method 

There are many methods for obtaining the minimum of the functional C knowing its gradient 
with respect to the controls. Adjoint methods involve the following steps: 

1. Start with a set cq of shape coefficients 

2. Enforce 8Ca = 0 and = 0 by Ending a U that satisfies the steady state Euler 
equations and boundary conditions 

3. Enforce 5Cu = 0 by finding a A that satisfies the costate equations and boundary 
conditions 

4. Calculate V a £, if it is 0 we have found the minimum, otherwise 

5. Update a with V„U, using a proper stepsize. 

6. Restart from 2 until V„U = 0. 

The need to repeat steps 2 and 3 above many times can become prohibitively expensive 
for geometrically complex configurations requiring computational power near the limits of 
present capabilities. 

Ta’asan proposed in [19] an efficient way of solving the optimization problem. The main 
observation is the following. The solution of the optimization problem lies on the intersection 
of the state, costate and design hypersurfeces (in the state, costate and design spaces). 
Gradient based methods (including adjoint formulations) can be viewed as marching along 
the intersection of the state and costate hypersurfaces. This is an expensive process since 
each step requires the solution of two PDEs. The idea of the pseudo-time method was to 
perform the marching on the design space while improving the distance to the other two. 
Compare Fig. 2 and 3. The cost of such an iteration per step is significantly smaller than 
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that of gradient based methods. Its convergence has been shown by Ta’asan in [19] to be 
independent of the number of design variables. 

The design equation for a wide class of problems, including the one considered here, is 
defined on the boundary only. Thus, it can be viewed as an extra boundary condition, 
and the design variables as the additional variables to solve for. In some cases the design 
equation can be solved for the design variables and a simple implementation of the above 
idea exists. In other cases the design equation, viewed as an equation for the design variables 
keeping the state and costate fixed, may be singular and a more involved implementation is 
required. This is the case for the problem considered here. In such cases it is necessary to 
solve for the design variables together with the state and costate variables in a small vicinity 
of the boundary S. 

Thus, at each step of computation on the entire held, the design equation is satisfied 
together with the boundary conditions for the state and costate equations. The solution on 
the entire held 0 affects the result of the optimization on S through the values of U and A 
on the auxiliary boundary T; see Fig. 4. 

The algorithm is as follows: 

1. Start with a tentative set of cq. 

2. March in time, on the entire held, the state equation a few steps. 

3. March in time, on the entire held, the costate equation a few steps. 

4. Solve in S the state equation with its boundary conditions, the costate equation with 
its boundary conditions and compute V a £. 

5. If V a £ = 0, restart from step 2, repeating steps 3 and 4 until the state and costate 
equations are converged on the entire held. Otherwise take a” +1 = a” + /(V a £) and 
go to 4. 

We took a'i +l = a” — aV a £, where a is a parameter. One could try to solve the problem 
on the boundary, in step 4 above, using a direct solver. The way we propose here has the 
advantage of being a simple modihcation of adjoint method, and therefore can be easily 
implemented. 

5 First optimization experiments 

We introduce a discrete grid dehned as (aq,y m ) = (x 0 + /Ax,y(0) + mAy) where Ax is 
constant and Ay is a constant fraction of the local height of the nozzle; see Fig. 5. 

The steady solution of the Euler equations is obtained with a time- dependent technique, 
in the frame of an explicit hnite volume code. The conservative variables U are computed 
at the cell centers, and the huxes F and G are evaluated at the cell interfaces using the 
approximate Riemann solver in [14]. Higher-order accuracy is achieved using an Essentially 
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Non-Oscillatory scheme [7]. The flow field values at cell interfaces, used as initial conditions 
for the Riemann problem, are reconstructed by means of a linear interpolation and using 
a minimod limiter. The amplitude of the integration step is chosen according to the CFL 
condition. 

The costate equations are discretized on the same grid presented above. Since they have 
no conservative form, the numerical solution is obtained using the finite-difference scheme 
proposed in [8]. 

The computations are performed on a 40 X 20 grid. Total pressure and total temperature 
at the inlet are taken to be unity and <r(0, y) = 0. At the outlet the static pressure depends 
on the test-case considered. For the lower wall ordinate y(0) we have 


y(&) 


( 0 if -0.5 < ;r < 0 

J2t= i a i U +1 { x — l) 2 if 0 < x < 1 

10 if 1 < a; < 1.5 


We try to recover the pressure distribution obtained with the Euler solver, corresponding 
to the set of shape coefficients a = (0,2, 2,0). This means that the functional is 0 at 
the minimum. The outlet pressure is such that the flow presents a relevant shock at re- 
compression. 

Figure 6 shows the functional values at each step of computation. Figure 7 shows the 
the convergence history of the state equation, computed to second-order accuracy, and the 
convergence history of the costate equation. Finally in Fig. 8, we present the starting 
pressure distribution and the one obtained at the end of the optimization procedure. 

The practicability of this approach depends on the rate of convergence to the minimum. In 
fact the state and costate equations converge to the steady solution with a less favorable rate 
compared to that of a simple analysis. It is easily seen that since the shape is changing, the 
flow field must change accordingly and so must the residuals. Figure 9 shows a comparison 
of the residuals for the state equations in the case of a simple analysis to the residuals in 
the optimization case. The convergence rate is 3 to 4 times slower in the optimization case. 
Considering the costate equations, the cost of the optimization procedure turns to be of the 
order of 10 analyses; using the first of the two algorithms presented in Section 4 the factor 
of proportion is 100 to 200 depending on the updating strategy used. The CPU time needed 
on a DEC 3000/500 is 18 minutes with the algorithm presented. For the first algorithm of 
Section 4, 6 hours of CPU time were needed. 

The present rate of convergence could be improved by changing the way of updating 
the grid. In fact, close to the minimum, the entire grid is perturbed to update only the 
boundary. We believe that, close to the minimum, the rate of convergence can be improved 
by updating only the boundary points of the grid. In fact, the small difference between 
the desired pressure and the one obtained is due to the fact that close to the minimum 
the convergence rate of the equations is reduced. Therefore the pressure p* is obtained 
asymptotically. 
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6 Optimal shape for compressible flows 

The following examples represent situations for which the optimal solution is not generated 
with the same algorithm we used to study the optimization problem. In the first case, we 
recover a pressure distribution known theoretically and compare the shape obtained with 
the theoretical one. 

The Ringleb flow (see [17]) is a two-dimensional steady compressible isentropic flow, where 
subsonic, transonic and supersonic regions are represented. It describes a 180°-turn of a 
compressible flow; all the exact values of the flow properties are given by simple formulas 
dependent on the stream function and on the Mach number. We consider the portion of 
the flow confined between two streamlines, which may be regarded as solid walls. The 
maximum Mach number on the bottom streamline is 1.6 and the minimum 0.8. On the top 
streamline the maximum Mach number is 0.8. The theoretical Mach number isocontours 
for such a flow are shown in Fig. 10. 

The Ringleb pressure distribution on the bottom wall is taken as the desired distribution 
that we want to achieve. 

The lower wall is described by the following parameterization: 

y(0) = r 0 + (n - r 0 ) J2 a i sin (* 77 W 57“) 

i=i “ 0 ° 

where r 0 is the distance, measured from the point of intersection of the lines from the inlet 
and the outlet, to the first point on the lower wall and r i is the distance to the last point. 
The angles 6 0 and 9\ are relative to the first and last point respectively, and are measured 
from the line from the inlet. 

In principle, if we try to recover the pressure distribution on the lower wall, the solution 
is out of the design space. We don’t have any a priori knowledge of the values that the ay 
will assume and how close to the desired pressure we can get. 

In Fig. 11 it is seen that no visible difference can be appreciated between the theoretical 
wall shape and the optimal shape found. The points representing the two solutions do not 
overlap since they are computed on two different grids. In this case the functional eq.(4) 
is 6.70 • 10 -5 after 500 iterations of steps 2. and 3. of the second algorithm proposed. See 
Fig. 12. The CPU time required for this case is about 20 minutes. 

In the second case considered, we are concerned with a convergent nozzle. The lower wall 
is represented by a parameterization similar to that above. The inlet Mach number is 2.2 
and the grid is 80 X 40. The starting configuration with pressure contourlines is shown in 
Fig. 13. 

A relevant shock is present in the flow held and our objective is to eliminate it by requiring 
a smooth compression at the lower wall. The smooth pressure distribution is not perfectly 
attained, as is seen in Fig. 14. Nevertheless the recompression appears to be smooth and 
the shock is eliminated from the how held (Fig. 15). These results are obtained after 200 
iterations in about 35 minutes of CPU. The functional is decreased from 3.75 • 10 -3 to 
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1.35 • 10 -4 . The computation has been pursued for 2000 iterations and the functional value 
remained unchanged. 

Finally, an experiment using 8 shape coefficients is performed. In Fig. 16 it is seen that 
the convergence rates of the state and costate equations are not affected. The functional 
minimum is therefore attained with the same number of iterations as in the case of 4 shape 
coefficients. 


7 Conclusions 

The pseudo-time method was applied to optimization problems governed by the Euler equa- 
tions in two dimensions. The problem of matching the pressure distribution to a desired 
one was considered, both for subsonic and supersonic flows. The rate of convergence to 
the minimum for the cases considered is 3 to 4 times slower compared to that of the ana- 
lysis problem. Results were obtained for Ringleb flow and a shockless recompression case. 
The algorithm could be implemented with no substantial changes to existing adjoint based 
codes. Numerical results indicate that the method converges at a rate which is independent 
of the number of design variables. The method offers a powerful and inexpensive tool for 
the study of non-intuitive configurations for aerodynamic design. 
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Figure 1: Model Problem. 
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Figure 2: Point A represents the desired optimum, point B l lip starting configuration. In 
the standard adjoint method, point A is reached by following a narrow path corresponding 
to the intersection of plane S and C . At each step along A — B, the state and costate 
equations are iterated many times. 
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Figure 3: In the new approach a new path, A — B , is taken lying on the plane T representing 
all solutions to the design equation and the boundary conditions of state and costate equa- 
tions. The computational cost of working on this plane is equivalent to solving a problem 
one space dimension less than the original problem. The solution of the state and costate 
equations is achieved only when point A is reached. 
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Figure 6: Logarithm of the functional versus the number of iterations on the entire held. 



Figure 7: Convergence history for state and costate equations. Logarithm of the residuals. 







Figure 8: Target pressure distribution and optimal one. Starting pressure: constant distri- 
bution at 0.7 reference value. 



Figure 9: Convergence history of state equations in a simple analysis compared to the 
convergence of the state equations in an optimization procedure. 
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Figure 12: Logarithm of the functional versus the number of iterations on the entire held. 



Figure 13: Convergent nozzle starting configuration: pressure isocontours. 
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Figure 16: Convergent nozzle. State and costate convergence history: 4 shape coefficients 
versus 8 shape coefficients. 
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